!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!> Representation of the distribution function.
!!
!! \n
!!
!! The distribution function is represented as a collection of \f$f\f$
!! and the electric field and potential. Notice that, since the
!! distribution function \f$f\f$ and the electric field (more
!! precisely, the potential) enter two different problems, namely the
!! Vlasov equation and the Poisson problem, we need to provide
!! separate vector space operations for each of them. Here, the vector
!! space operations for the distribution function are attached to \c
!! t_f_state, while a separate type is introduced for the electric
!! potential and the associated Lagrange multiplier.
!!
!! Summarizing, the following types are provided:
!! <ul>
!!  <li> \c t_f_state: prognostic variable \f$f\f$, extended from \c
!!  c_stv and thus with the corresponding vector space operations
!!  <li> \c t_e_field: diagnostic variables related with the
!!  electromagnetic field \f$\Phi\f$, \f$\underline{E}\f$,
!!  \f$\tilde{\underline{E}}\f$ and \f$\rho\f$, extended from \c c_ods
!!  <ul>
!!   <li> \c t_p_state: used to represent \f$\Phi\f$ providing the
!!   vector space operations required by the linear solvers, since
!!   this is the solution of the Poisson problem (see also \c
!!   mod_poisson_pb).
!!  </ul>
!! </ul>
!<----------------------------------------------------------------------
module mod_f_state

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_octave_io, only: &
   mod_octave_io_initialized, &
   write_octave

 use mod_state_vars, only: &
   mod_state_vars_initialized, &
   c_stv

 use mod_time_integrators, only: &
   mod_time_integrators_initialized, &
   c_ode, c_ods

 use mod_tps_phs_grid, only: &
   mod_tps_phs_grid_initialized, &
   t_tps_phs_grid

 use mod_tps_base, only: &
   mod_tps_base_initialized, &
   t_tps_base

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_f_state_constructor, &
   mod_f_state_destructor,  &
   mod_f_state_initialized, &
   t_p_state, t_f_state, t_e_field, &
   new_f_state, new_e_field, clear, &
   write_octave

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !> Electric potential
 !!
 !! This type extends \c c_stv because it is used in the linear
 !! solvers.
 type, extends(c_stv) :: t_p_state
  !> The degrees of freedom are stored in a 1D array because this is
  !! more convenient for the solution of the Poisson problem. The last
  !! element of the array is the Lagrange multiplier.
  real(wp), allocatable :: pl(:)
 contains
  procedure, pass(x) :: incr => p_incr
  procedure, pass(x) :: tims => p_tims
  procedure, pass(x) :: inlt => p_inlt
  procedure, pass(z) :: alt  => p_alt
  procedure, pass(z) :: mlt  => p_mlt
  procedure, pass(z) :: copy => p_copy
  procedure, pass(x) :: scal => p_scal
  procedure, pass(x) :: show => p_show
  ! To be removed when the corresponding compiler bug are resolved
  procedure, pass(x) :: source       => p_source
  procedure, pass(x) :: source_vect  => p_source_vect
 end type t_p_state

 !> State information for the distribution function
 !!
 !! Some pointers to the finite element space description are also
 !! introduced; such description is shared by all the state
 !! variables.
 type, extends(c_stv) :: t_f_state
  type(t_tps_phs_grid), pointer :: grid    => null()
  type(t_tps_base),     pointer :: base(:) => null()
  !> The distribution function is organized on three indexes: the
  !! first one identifies the local x degrees of freedom, the second
  !! one the local v degrees of freedom, and the third one the
  !! elements.
  real(wp), allocatable :: f(:,:,:)
 contains
  procedure, pass(x) :: incr
  procedure, pass(x) :: tims
  !> This subroutine is overridden both because of efficiency and
  !! because many compilers leak memory when working with polymorphic
  !! temporaries
  procedure, pass(x) :: inlt => f_inlt
  procedure, pass(z) :: alt  => f_alt
  procedure, pass(z) :: mlt  => f_mlt
  procedure, pass(z) :: copy => f_copy
  procedure, pass(x) :: scal => f_scal
  procedure, pass(x) :: show => f_show
  procedure, pass(x) :: source
  procedure, pass(x) :: source_vect
  !> The next subroutine is useful for debug
 end type t_f_state

 !> Diagnostics: electric field and potential
 type, extends(c_ods) :: t_e_field
  type(t_tps_phs_grid), pointer :: grid    => null()
  type(t_tps_base),     pointer :: base(:) => null()
  !> The electric field is organized on three indexes: the spatial
  !! component, the local degree of freedom and the space element
  real(wp), allocatable :: e(:,:,:)
  !> The projected electric field, expressed in terms of the scalar
  !! basis in space variable. An additional index (the second one) is
  !! introduced because various projections are used, depending on the
  !! sign of the velocity components. In particular, each component of
  !! the projected electric field comes in two <em>flavours</em>.
  real(wp), allocatable :: et(:,:,:,:)
  !> Density (same indexes as \c p)
  real(wp), allocatable :: r(:,:)
  !> Electric potential
  type(t_p_state) :: p
 ! private members
 end type t_e_field

 interface clear
   module procedure clear_f_state, clear_e_field
 end interface

 interface write_octave
   module procedure write_p_state, write_f_state, write_e_field
 end interface write_octave

! Module variables

 ! public members
 logical, protected ::               &
   mod_f_state_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_f_state'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_f_state_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_kinds_initialized.eqv..false.) .or. &
    (mod_messages_initialized.eqv..false.) .or. &
   (mod_octave_io_initialized.eqv..false.) .or. &
  (mod_state_vars_initialized.eqv..false.) .or. &
(mod_time_integrators_initialized.eqv..false.) .or. &
(mod_tps_phs_grid_initialized.eqv..false.) .or. &
    (mod_tps_base_initialized.eqv..false.) ) then  
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_f_state_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_f_state_initialized = .true.
 end subroutine mod_f_state_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_f_state_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_f_state_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_f_state_initialized = .false.
 end subroutine mod_f_state_destructor

!-----------------------------------------------------------------------
 
 subroutine new_f_state(obj,grid,base)
  type(t_tps_phs_grid), intent(in), target :: grid
  type(t_tps_base),     intent(in), target :: base(2)
  type(t_f_state),      intent(out) :: obj

   obj%grid => grid
   obj%base => base
   allocate( obj%f(base(1)%pkd,  base(2)%pkd ,grid%ne   ) )

 end subroutine new_f_state
 
!-----------------------------------------------------------------------
 
 pure subroutine clear_f_state(obj)
  type(t_f_state), intent(out) :: obj

   nullify( obj%grid )
   nullify( obj%base )
   ! allocatable components are implicitly deallocated
 
 end subroutine clear_f_state
 
!-----------------------------------------------------------------------

 subroutine new_e_field(obj,grid,base)
  type(t_tps_phs_grid), intent(in), target :: grid
  type(t_tps_base),     intent(in), target :: base(2)
  type(t_e_field),      intent(out) :: obj

   obj%grid => grid
   obj%base => base
   allocate( obj%e(   grid%d,  base(1)%pkdx,grid%gx%ne) )
   allocate( obj%et(  grid%d,2,base(1)%pkd ,grid%gx%ne) )
   allocate( obj%p%pl(         base(1)%pkd*grid%gx%ne+1))
   allocate( obj%r(            base(1)%pkd ,grid%gx%ne) )

 end subroutine new_e_field
 
!-----------------------------------------------------------------------
 
 pure subroutine clear_e_field(obj)
  type(t_e_field), intent(out) :: obj

   nullify( obj%grid )
   nullify( obj%base )
   ! allocatable components are implicitly deallocated
 
 end subroutine clear_e_field
 
!-----------------------------------------------------------------------

 subroutine p_incr(x,y)
  class(c_stv),     intent(in)    :: y
  class(t_p_state), intent(inout) :: x

   select type(y); type is(t_p_state)

   x%pl = x%pl + y%pl

   end select
 end subroutine p_incr

!-----------------------------------------------------------------------

 subroutine incr(x,y)
  class(c_stv),     intent(in)    :: y
  class(t_f_state), intent(inout) :: x

   select type(y); type is(t_f_state)

   x%f = x%f + y%f

   end select
 end subroutine incr

!-----------------------------------------------------------------------

 subroutine p_tims(x,r)
  real(wp),         intent(in)    :: r
  class(t_p_state), intent(inout) :: x

   x%pl = r*x%pl

 end subroutine p_tims

!-----------------------------------------------------------------------

 subroutine tims(x,r)
  real(wp),         intent(in)    :: r
  class(t_f_state), intent(inout) :: x

   x%f = r*x%f

 end subroutine tims

!-----------------------------------------------------------------------

 subroutine p_inlt(x,r,y)
  real(wp),     intent(in) :: r
  class(c_stv), intent(in) :: y
  class(t_p_state), intent(inout) :: x

   select type(y); type is(t_p_state)

   x%pl = x%pl + r*y%pl

   end select
 end subroutine p_inlt

!-----------------------------------------------------------------------

 subroutine f_inlt(x,r,y)
  real(wp),     intent(in) :: r
  class(c_stv), intent(in) :: y
  class(t_f_state), intent(inout) :: x

   select type(y); type is(t_f_state)

   x%f = x%f + r*y%f

   end select
 end subroutine f_inlt

!-----------------------------------------------------------------------

 subroutine p_alt(z,x,r,y)
  real(wp),     intent(in) :: r
  class(c_stv), intent(in) :: x, y
  class(t_p_state), intent(inout) :: z

   select type(x); type is(t_p_state)
    select type(y); type is(t_p_state)

   z%pl = x%pl + r*y%pl

    end select
   end select
 end subroutine p_alt

!-----------------------------------------------------------------------

 !> Overloaded for efficiency
 subroutine f_alt(z,x,r,y)
  real(wp),     intent(in) :: r
  class(c_stv), intent(in) :: x, y
  class(t_f_state), intent(inout) :: z

   select type(x); type is(t_f_state)
    select type(y); type is(t_f_state)

   z%f = x%f + r*y%f

    end select
   end select
 end subroutine f_alt

!-----------------------------------------------------------------------

 subroutine p_mlt(z,r,x)
  real(wp),     intent(in) :: r
  class(c_stv), intent(in) :: x
  class(t_p_state), intent(inout) :: z

   select type(x); type is(t_p_state)

   z%pl = r*x%pl

   end select
 end subroutine p_mlt

!-----------------------------------------------------------------------

 !> Overloaded for efficiency
 subroutine f_mlt(z,r,x)
  real(wp),     intent(in) :: r
  class(c_stv), intent(in) :: x
  class(t_f_state), intent(inout) :: z

   select type(x); type is(t_f_state)

   z%f = r*x%f

   end select
 end subroutine f_mlt

!-----------------------------------------------------------------------

 subroutine p_copy(z,x)
  class(c_stv),     intent(in)    :: x
  class(t_p_state), intent(inout) :: z

   select type(x); type is(t_p_state)

   z%pl = x%pl

   end select
 end subroutine p_copy

!-----------------------------------------------------------------------

 subroutine f_copy(z,x)
  class(c_stv),     intent(in)    :: x
  class(t_f_state), intent(inout) :: z

   select type(x); type is(t_f_state)

   z%f  = x%f

   end select
 end subroutine f_copy

!-----------------------------------------------------------------------

 function p_scal(x,y) result(s)
  class(t_p_state), intent(in) :: x
  class(c_stv), intent(in) :: y
  real(wp) :: s
   
   select type(y); type is(t_p_state)
    s = dot_product(x%pl,y%pl)
   end select

 end function p_scal

!-----------------------------------------------------------------------

 function f_scal(x,y) result(s)
  class(t_f_state), intent(in) :: x
  class(c_stv), intent(in) :: y
  real(wp) :: s
   
   select type(y); type is(t_f_state)
    s = sum(x%f*y%f)
   end select

 end function f_scal

!-----------------------------------------------------------------------

 subroutine p_source(y,x)
  class(t_p_state), intent(in) :: x
  class(c_stv), allocatable, intent(out) :: y

   allocate(t_p_state::y)
 
   select type(y); type is(t_p_state)

   allocate(y%pl(size(x%pl)))

   end select
 end subroutine p_source
 subroutine p_source_vect(y,x,m)
  integer, intent(in) :: m
  class(t_p_state), intent(in)  :: x
  class(c_stv), allocatable, intent(out) :: y(:)

  integer :: i

   allocate(t_p_state::y(m))

   select type(y); type is(t_p_state)

   do i=1,m
     allocate(y(i)%pl(size(x%pl)))
   enddo

   end select
 end subroutine p_source_vect

!-----------------------------------------------------------------------

 subroutine p_show(x)
  class(t_p_state), intent(in) :: x
   write(*,'(a)') '------- Show a t_p_state object -------'
   write(*,'(*(e12.3))') x%pl
   write(*,'(a)') '---------------------------------------'
 end subroutine p_show

!-----------------------------------------------------------------------

 subroutine f_show(x)
  class(t_f_state), intent(in) :: x
   write(*,'(a)') '------- Show a t_f_state object -------'
   write(*,'(*(e12.3))') maxval(abs(x%f))
   write(*,'(a)') '---------------------------------------'
 end subroutine f_show

!-----------------------------------------------------------------------

 subroutine source(y,x)
  class(t_f_state), intent(in) :: x
  class(c_stv), allocatable, intent(out) :: y

   allocate(t_f_state::y)
 
   select type(y); type is(t_f_state)

   call new_f_state( y , x%grid , x%base )

   end select
 end subroutine source
 subroutine source_vect(y,x,m)
  integer, intent(in) :: m
  class(t_f_state), intent(in)  :: x
  class(c_stv), allocatable, intent(out) :: y(:)

  integer :: i

   allocate(t_f_state::y(m))

   select type(y); type is(t_f_state)

   do i=1,m
     call new_f_state( y(i) , x%grid , x%base )
   enddo

   end select
 end subroutine source_vect

!-----------------------------------------------------------------------

 subroutine write_p_state(p,var_name,fu)
  integer, intent(in) :: fu
  type(t_p_state), intent(in) :: p
  character(len=*), intent(in) :: var_name
 
  character(len=*), parameter :: &
    this_sub_name = 'write_p_state'

   write(fu,'(a,a)') '# name: ',var_name
   write(fu,'(a)')   '# type: struct'
   write(fu,'(a)')   '# length: 1' ! number of fields

   ! field 01 : pl
   write(fu,'(a)')   '# name: pl'
   write(fu,'(a)')   '# type: cell'
   write(fu,'(a)')   '# rows: 1'
   write(fu,'(a)')   '# columns: 1'
   call write_octave(p%pl,'c','<cell-element>',fu)

 end subroutine write_p_state

!-----------------------------------------------------------------------

 subroutine write_f_state(f,var_name,fu)
  integer, intent(in) :: fu
  type(t_f_state), intent(in) :: f
  character(len=*), intent(in) :: var_name
 
  character(len=*), parameter :: &
    this_sub_name = 'write_f_state'

   write(fu,'(a,a)') '# name: ',var_name
   write(fu,'(a)')   '# type: struct'
   write(fu,'(a)')   '# length: 2' ! number of fields

   ! field 01 : f
   write(fu,'(a)')   '# name: f'
   write(fu,'(a)')   '# type: cell'
   write(fu,'(a)')   '# rows: 1'
   write(fu,'(a)')   '# columns: 1'
   call write_octave(f%f,'<cell-element>',fu)

 end subroutine write_f_state

!-----------------------------------------------------------------------

 subroutine write_e_field(ef,var_name,fu)
  integer, intent(in) :: fu
  type(t_e_field), intent(in) :: ef
  character(len=*), intent(in) :: var_name
 
  character(len=*), parameter :: &
    this_sub_name = 'write_e_field'

   write(fu,'(a,a)') '# name: ',var_name
   write(fu,'(a)')   '# type: struct'
   write(fu,'(a)')   '# length: 4' ! number of fields

   ! field 01 : e
   write(fu,'(a)')   '# name: e'
   write(fu,'(a)')   '# type: cell'
   write(fu,'(a)')   '# rows: 1'
   write(fu,'(a)')   '# columns: 1'
   call write_octave(ef%e,'<cell-element>',fu)

   ! field 02 : et
   write(fu,'(a)')   '# name: et'
   write(fu,'(a)')   '# type: cell'
   write(fu,'(a)')   '# rows: 1'
   write(fu,'(a)')   '# columns: 1'
   call write_octave(ef%et,'<cell-element>',fu)

   ! field 03 : p
   write(fu,'(a)')   '# name: p'
   write(fu,'(a)')   '# type: cell'
   write(fu,'(a)')   '# rows: 1'
   write(fu,'(a)')   '# columns: 1'
   call write_octave(ef%p,'<cell-element>',fu)

   ! field 04 : r
   write(fu,'(a)')   '# name: r'
   write(fu,'(a)')   '# type: cell'
   write(fu,'(a)')   '# rows: 1'
   write(fu,'(a)')   '# columns: 1'
   call write_octave(ef%r,'<cell-element>',fu)

 end subroutine write_e_field

!-----------------------------------------------------------------------

end module mod_f_state

